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We study the effect of strong disorder in a 3-dimensional topological insulators with time-reversal 
symmetry and broken inversion symmetry. Firstly, using level statistics analysis, we demonstrate 
the persistence of delocalized bulk states even at large disorder. The delocalized spectrum is seen 
to display the levitation and pair annihilation effect, indicating that the delocalized states continue 
to carry the Z2 invariant after the onset of disorder. Secondly, the Z2 invariant is computed via 
twisted boundary conditions using a novel and efficient numerical algorithm. We demonstrate that 
the Z2 invariant remains quantized and non-fluctuating even after the spectral gap becomes filled 
with dense localized states. In fact, our results indicate that the Z2 invariant remains quantized 
until the mobility gap closes or until the Fermi level touches the mobility edges. Based on such data, 
we compute the phase diagram of the Bi2Se3 topological material as function of disorder strength 
and position of the Fermi level. 
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Topological insulators represent a new class of mate- 
rials where the topology of the bulk electronic structure 
induces highly non-trivial effects. One such effect is 
the emergence of metallic states along the edges of pla- 
nar (2D) topological insulating structures or on the sur- 
face of 3-dimensional (3D) topological insulators. These 
materials became a reality after a topological insulator 
from the class of Quantum Spin-Hall insulators and one 
from the class of strong 3D topological insulators with 
time reversal symmetry have been theoretically predicted 
and then engineered and characterized in labor at ories.^^^ 
Since then, many additional topological materials have 
been discovered (the reader can find a survey of the field 
in Refs. 9-12). 

For topological insulators with time-reversal symme- 
try, it has been argued that the edge/surface states main- 
tain their metallic character even in the presence of weak 
disorder, ^'^ due to the cancelation of the backscattering 
amplitudes. This robustness against disorder can be the 
key to many technological applications,^"^ ^ and because 
of that, a great deal of effort has been dedicated to under- 
standing the behavior of the topological materials in the 
presence of disorder. ^^^^ One important question, which 
is still opened for 3D topological insulators, is if the ro- 
bustness against disorder extends into the strong disorder 
regime, particularly into the regime where the insulating 
gap is filled with dense localized spectrum. The theo- 
retical argument based on the time-reversal symmetry is 
perturbative and therefore it breaks down in this regime. 
As such, one must seek a new argument that combines 
topology (an index theorem) and symmetry and this has 
become an extremely active area of research. While 
searching for such argument, it became apparent to us 
that a numerical exploration of the matter will be of great 
help. Notable efforts paralleling ours are the numerical 
computations of a newly defined Bott index in Ref. 28, 
carried out for a disordered model of a 3D strong topo- 



logical insulator. The explicit equivalence between the 
Bott index and the strong Z2 invariant remains to be es- 
tablished in the strong disorder regime considered here. 

Our discussion will be restricted to 3D topological insu- 
lators. In the presence of time reversal symmetry, the in- 
sulators follow a Z2 topological classification. The strong 
Z2 invariant that renders an insulator as either trivial or 
topological was computed in various ways, but in general 
the computations were quite demanding because they 
had to be carried out with special smooth gauges. The 
difficulty introduced by this requirement has been docu- 
mented in our previous work.^^ For example, the origi- 
nal expressions of the Z2 invariant, ^'^'^^'^^ require special 
smooth gauges and were computed only for analytically 
solvable band models. These expressions have been refor- 
mulated in an almost gauge invariant fashion by Fukui 
and his collaborators. ^^'^^ The method still requires a 
time-reversal adapted gauge at the boundary of half of 
the Brillouin zone, but nevertheless it became the method 
of choice when computing the Z2 invariant. ^^'^^^^ Still, 
an application of the method to the disordered case ex- 
ists only in 2D.^^ The Chern-Simons integral of the quan- 
tized magneto-electric polarization also requires a glob- 
ally smooth gauge. The difficulties introduced by this 
requirement were highlighted in Ref. 40, where the best 
effort to evaluate the Chern-Simons integral only led to a 
value of 0.3, for a disorder- free topological case where the 
result should have been quantized to 1. To date, there 
is no successful direct evaluation of this Chern-Simons 
integral for clean tight-binding models. 

The issue was recently reconsidered and gauge- 
independent formulations of the weak and strong Z2 in- 
variants are now available. ^^'^^'^^^^ Here we will follow 
Ref. 29 and we will argue here that this new formu- 
lations bring certain numerical advantages which open 
the possibility of directly computing the Z2 invariants 
for systems with extremely large unit cells, particularly 
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for disordered samples (as opposed to indirectly inferring 
the Z2 invariants from other type of calculations such as 
transport simulations of the surface states). We present a 
numerical analysis of the strong invariant for a system 
without inversion symmetry, computed in the weak and 
strong disorder regimes via the twisted boundary condi- 
tions technique combined with the new formulation of the 
invariant. The use of the twisted boundary conditions 
was advocated by Kane and Mele in their original dis- 
cussion of the 2D Z2 invariant as an effective procedure 
for tackling the effect of disorder and electron- electron 
interaction.^ Numerically, this method is equivalent to 
computing the Z2 invariant for a periodic system with 
a very large unit cell, leading to thousands of occupied 
energy bands. Finding smooth special gauges for such 
complex band structures is prohibitively difficult, which 
is why the Z^ invariant has never been computed be- 
fore for disordered 3D topological insulators (the parity 
analysis was accomplished in Ref. 23 for a system and 
disorder with inversion symmetry). 

Working with a tight-binding model appropriate for 
the 3D topological material Bi2Se3, we provide com- 
pelling evidence that the strong Z^ invariant remains 
well defined and quantized even after the insulating gap 
becomes ffiled with dense localized spectrum. In fact, 
our various mappings of the Z2 invariant indicate that 
the quantization holds as long as the Fermi level remains 
in the mobility gap. Furthermore, we use level statis- 
tics analysis to map the localized or extended character 
of the energy spectrum for a wide range of disordered 
strengths. We provide compelling evidence that there are 
bulk metallic states that persist even at very large disor- 
der and we derive the phase diagram of the 3D model as 
function of Fermi level and disorder strength. The phase 
diagram consists of the strong topological phase, which 
is completely surrounded by a metallic phase, which is 
again surrounded by the trivial insulating phase. Com- 
putations of the Z2 invariant along paths that cross from 
the topological into the trivial insulating phase reveal 
strict quantization of the invariant to ±1 values in the 
topological/trivial insulating phases, respectively, and 
strong fluctuations between ±1 inside the metallic phase. 

The motivation behind the present study was two-fold. 
Firstly, there is no theory of the Z2 invariant for aperi- 
odic systems (except for the trivial case when the Fermi 
level is in a spectral gap, i.e. a region that is void of 
any energy spectrum). For example, for the Chern and 
spin-Chern invariants we have theories that provide ex- 
plicit and specific conditions, which can be written in 
one line, that tells us when these invariants take quan- 
tized values even if the Fermi level is not in a spectral 
gap.^^'^^ Furthermore, we have real-space formulations 
of these two invariant s,^^'^^ which allows one to compute 
them in "one shot" without involving twisted boundary 
conditions. Nothing like that exists for the strong Z^ in- 
variant, despite some sustained efforts. This makes one 
to question that such a theory will ever be achieved, and 
in fact to question that the strong Z2 invariant does in- 



deed remain quantized once the spectral gap is closed. 
Our numerical study provides the first direct evidence 
that the strong Z^ invariant behaves similarly with these 
other two invariants, which can be a strong motivation 
for people to continue searching for a theory of the strong 
Z2 invariant for aperiodic systems. 

Secondly, it is known that disorder can strongly deform 
the phase boundary of the topological state. -^^'^^'^^'^^ As 
such, it is highly desired to devise quantitative methods 
that can accurately pinpoint the extent of the topologi- 
cal phases. Previously, the strong topological phase was 
identified by probing the metallic character of the surface 
states via transport calculations in a long bar geometry. 
For the special case when the system and the disorder 
have inversion symmetry, the topological phase was iden- 
tified using the parities of the states. Ideally, however, 
it will be to directly map the strong Z^ invariant and 
our study demonstrates that this is indeed possible for 
3D materials. 

Lastly, we want to state explicitly that our numerical 
simulations probe un-charted territories. As discussed in 
the next sections, our method, and for that matter all the 
established methods, are easily seen to produce quantized 
and non- fluctuating Z2 values, if a spectral gap between 
the occupied and non-occupied levels remains open at all 
times while twisting of the boundary conditions. How- 
ever, at strong disorder, the spectral gap not only closes 
but the levels can change their ordering when twisting 
the boundary conditions. Thus, levels that once were 
occupied become un-occupied and vice-versa. Moreover, 
different disorder configurations can no longer be con- 
nected adiabatically. While our numerical procedure can 
still be applied in these situations, the available theoret- 
ical arguments can no longer assure us that the output 
remains the same from one disorder configuration to an- 
other (this also applies to the Bott index of Refs. 27 and 
28). Still, if the states near the Fermi level are local- 
ized, one expects the value of the invariant to remain 
un-affected by this phenomenon. This is exactly what 
we are trying to verify in this work. 



I. THE TWISTED BOUNDARY CONDITIONS 

In the first half of the paper we will carry the discussion 
at a general level. We consider a generic 3D quantum 
lattice model with many quantum states ^ per site n. 
The Hilbert space 1-i is spanned by |n, ^) and the periodic 
Hamiltonian is given by: 

Ho= Yl \n,0hf{n+p,^'\=J2\r^)hp{n+P\,{^) 

n,p,^,^' n,p 

where p runs over first, second, etc., neighbors of the ori- 
gin, \n) denotes the one column matrix with the entries 
|n,^), (n| represents its dual, and hp is the matrix of 
elements . We are interested in the properties of a 
disordered Hamiltonian = Hq-\-V^, and to be specific 
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we consider an onsite random potential: 



V^ = wY,^r^^W.O{n.^\ = wY, {n\ (2) 



n 



where cj^i^ are randomly independent amplitudes uni- 
formly distributed in the interval [— ^, 

The twisted boundary conditions method for disor- 
dered systems consists basically in considering a large 
super-cell 5, containing many unit cells of the clean sys- 
tem, specifically the states |n, < ^1,722,^3 < A/" — 1, 
and a random potential is placed inside this super-cell. 
The super-cell is then periodically repeated in space. We 
will continue to use to denote the resulting approx- 
imate Hamiltonian. Since we are dealing with a peri- 
odic system, we can construct a dual /c-space representa- 
tion via the Bloch transformation. This transformation is 
given by the isometry U from the Hilbert space 1-L of the 
infinitely repeated system into a continuum direct sum 
of copies of the supercell's Hilbert space Ho spanned by 
|n, < ni, 77-2, ns < A/" — 1 (T = 3D torus): 



U'.n 

U\n^mM,Cl 



e Ho, 

fcer 

e 

ker 



(27r)3/2 



(3) 



for any n in S and arbitrary m in Z^. Note that any 
point from can be uniquely written as n + mj\f. The 
inverse of the isometry is: 



k'er 



(27r) 



1 ^ e^^-^\n^m^,0- 



(4) 



Note that we wrote the transformation so that all the 
states inside the supercell get the same phase factor, in 
which case the 27r-periodicity in the /c-variables is auto- 
matically satisfied. Under this transformation we have: 



(5) 



where the Bloch Hamiltonians H^^{k) : Ho ^ Ho are 
defined by: 



nes p 

+W E \n)v^{n\ 
nes 



(6) 



with 



Kpik) = e »=i hp 



(7) 



These are precisely the twisted boundary conditions since 
the phase factor above occurs only for the lattice points 
n at the boundary. At this point we obtained a family of 
Bloch Hamiltonians indexed by a point on the 3D torus. 
The same construction can be achieved by wrapping the 
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FIG. 1. Example of a time-reversal invariant path in the 
Brillouin torus and its discretization. 



supercell S into a 3D torus and by threading magnetic 
fluxes through the 2D sections of this torus. The effect 
of such magnetic fluxes is captured by the same twisted 
boundary conditions. 

One should note that the twisted boundary conditions 
method is defined for a finite volume, whereas one is 
really interested in the infinite bulk samples. Strictly 
speaking, one has to take the volume of the supercell to 
infinity and carefully investigate the stability of the re- 
sults. In practice, of course, we will have to stop the limit 
at some point. 



II. THE Z2 INVARIANTS 

We start our discussion of the Z2 invariant with a brief 
review of the formulation given in Ref. 29. The connec- 
tion between this formulation and the previously exist- 
ing ones has been exhaustively discussed in this reference 
and will not be addressed here. Instead, we will give a 
detailed discussion of the numerical advantages brought 
in by this new method. 

Briefiy, the method goes as follows. Let Pfc denote 
the projector on to the states of H^{k) below the Fermi 
level Ep^ and let denote the time-reversal operation 
and assume the time- reversal invariance: 



(8) 



One considers a closed, time-reversal invariant path (i.e. 
a path which is mapped into itself by 6) on the Brillouin 
torus, 



[-7r,7r] 3 k ^ k{k), 



(9) 



parametrized by the variable k (we chose the notation 
on purpose because in practice this variable will often 
be kz for example). Then one integrates the differential 
equation (with the initial condition Uk',k' = Pk{k'))'' 



(10) 



We will use the simplified notation: Pk{k) = Pk- The 
result of the integration gives the unitary time evolution 



operator Uk^k' corresponding to the process of adiabat- 
ically changing the /c-vector along the chosen path in 
the Brihouin torus. It is assumed that the path starts 
{k = — tt) and closes {k = tt) at a time-reversal invariant 
fc-point. Necessarily, the path will cross another time- 
reversal invariant point at midway k = (see Fig. 1). 
Next, one considers arbitrary bases {e^} and {e^} for 
the occupied spaces at the time-reversal invariant points 
A: = and /c = tt, respectively, and one defines the fol- 
lowing matrices: 

0% = {el\e\e^p), (11) 

These matrices satisfy the following fundamental 
relation: 

PfR}-Met{?7}Pf{^o} ^ ^^2^ 

The left hand side of Eq. 12 will be called a pseudo 
Z2 invariant for the following reasons. The left hand 
side is gauge independent. Indeed, given the transfor- 
mation properties of the Pfaffians and determinants un- 
der the conjugation with unitary matrices, one can easily 
see that the numerator is independent of the bases {e^} 
and {ej}.^^ At the denominator, inside the square root, 
U-j^^-TT maps the k = ±7r occupied space into itself, so at 
a change of {ej} basis we have U^^-t^ EUt^^-t^E~'^ ^ 
with E a unitary matrix, so the determinant remains un- 
changed. However, the sign in Eq. 12 depends on which 
branch of the square root is used, but once a choice is 
made the value of the left-hand side cannot be changed 
by smooth deformations of the Bloch Hamiltonians that 
keep the insulating gap opened. 

In practice, the adiabatic evolution operators is com- 
puted by discretizing the paths and taking the product 
of projectors onto the occupied spaces at these discrete 
/c-points. Since the path is time-reversal invariant, we 
can choose the discretization points so as feg, A^i, . . . , fe^ 
discretizes the path from A: = to = tt, while —k^ 
—kn-i, ...,ko discretizes the path from k = — tt to A: = 0. 
In this case: 

U^,-^= lim Pfc^Pfc_, ...Pfc,...P_fc_,P_fc^. (13) 

In practice however, we have to stop limit at some n = n 
and work with an approximation: 

f/.,-. = Pk^Pk^_, . . . Pfeo • • • P-k^-iP-k^ , (14) 
and similar for U^^^o- 

U.fi = Pk^Pk^.^...Pko- (15) 

We are going to show in the following that the quantiza- 
tion in Eq. 12 remains exact even for finite n's. 
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Indeed, using the elementary fact that 6Pk.6~^ = 
P_fe^. , we have: 

U^,-n = PkA-i • • • PkoOPko • • • Pk^^.PkJ-^ (16) 

Inserting the identity operator J2a \^a){^a\ appro- 
priate places, we obtain: 

{el\U^,-^\eJ) = (e-|Pfe,Pfe,_, ...PfeJeO)(eO|0|eO) 

x(eO|Pfe„...Pfe„_,PfeJe|)(e||^-i|ep 

= {el\PkA.,...Pko\e1){eo)sy 

x{eJ\PkA-.---PkoK){0^')ip. 

(17) 

Summation over repeating indices was assumed above. 
At this step, the conclusion is: 

U^_^ = U§oU^O-\ (18) 

Taking the determinant and using the elementary prop- 
erties of the determinants and pfaffians we obtain: 

det{[/^,_^} = [Pf{4}"' det{f7}Pf{^o}]', (19) 

which is precisely Eq. 12. 

The significance of the above conclusion for the numer- 
ical calculations is that it allows us to use relatively small 
number of discretization points when evaluating Eq. 12. 
One question that could be asked is if the result of such 
calculation, while indeed quantized, it really equals the 
result in the n ^ 00 limit? To answer this question, 
we imagine a calculation with a dense number of dis- 
cretization points and then adiabatically collapsing pairs 
of adjacent discretization points into a single discretiza- 
tion point. In this way we can adiabatically transform 
the original computation into a computation with half 
the number of discretized points. Repeating the same 
action we can adiabatically reduce the number of dis- 
cretization points even further, by 4, 8 and so on. Since 
Eq. 12 is quantized, it cannot change its value during 
such adiabatic deformations, if all the quantities remain 
well defined. So what can go wrong? In the n ^ 00 limit, 
Uj^^-n is a true unitary operator so its determinant is a 
complex number on the unit circle. For finite n, Ut^^-tt is 
no longer unitary and its determinant moves inside the 
unit circle. As the number of discretization points is re- 
duced, the determinant moves closer to the origin so there 
is the possibility that det{UT^^-^} actually becomes equal 
to zero. At such instance, the calculation breaks down 
and the quantized value of Eq. 12 can change. So the 
conclusion is that Eq. 12 can be indeed evaluated using a 
relatively small number of discretization points, as long 
as one makes sure that det{f/7^^_7^} does not touches the 
origin. In practice we choose the number of discretized 
points so that | det{C/7P^_7p}| ^ 0.5, which reduces the 
number of required A:-points by an order of magnitude 
in our calculations, when compare with the case when 
|det{C/,,_,}|^0.9. 
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Eq. 12 is fundamental for the formulation of the Z2 in- 
variant but it cannot define an invariant by itself. That 
is because we don't have a canonical way to choose the 
branch of the square root at the denominator of Eq. 12. 
However, the important observation is that if one con- 
siders a pair of paths, then there is a canonical way to 
choose the same branch of the square root and genuine 
Z2 invariants can be defined. This has been detailed in 
Ref. 29. The following lines explain how the procedures 
were explicitly implemented in our calculations. 

For a 3D system, we consider 4 independent time- 
reversal invariant paths. If Vk^.ky denotes the path along 
kz direction that intersects the plane /C;^ = at {kx^ky), 
then we choose the following 4 paths: 

Po,o: (0,0,-^)^(0,0,^) 

Vo,n ' (0,7r, -tt) ^ (0,7r,7r) , . 

V^,o'- (^,0,-^)^(^,0,^) ^^^^ 

7^7r,7r : (7r,7r, -tt) ^ (7r,7r,7r), 

We interpolated between the paths 7^0,0 and Vo,7r using 
the process: 

[0,^] 3 ky^Vo^ky (21) 
By computing the adiabatic evolution 

U(0,ky-7T)^{0,ky,7r) (22) 

for the path Vo,ky^ we continuously interpolate between 
the determinants 

det{/7(o,0,-7r)^(0,0,7r)} ^ det{[/(o,7r,-7r)^(0,7r,7r)}: (23) 

using the process: 

[0,7r] 3ky^ Ciet{[/(0,fc„-7r)^(0,/c„7r)}. (24) 

This allows us to monitor how the determinant moves on 
the Riemann surface of the square root function, and to 
determine the location of det{[/(o,7r,-7r)^(o,7r,7r)} relative 
to the location of det{[/(o,o,-7r)^(o,o,7r)} the Riemann 
surface. If det{[/(o,/c^,-7r)^(o,/c^,7r)} crosses the semi-axis 
(— oc, 0) an odd number of times, then the determinants 
are located on opposite Riemann sheets and we have to 
use different branches of the square root, that is, we will 
have to use zby^ for one determinant and =F\/^ for the 
other determinant in Eq. 23 when we evaluate the de- 
nominator of Eq. 12. If det{[/(o,/c^,-7r)^(o,/cy,7r)} crosses 
the semi- axis (—00, 0) an even number of times, then the 
determinants are located on the same Riemann sheet and 
we have to use =by^ for one determinant and same dzy^ 
for the other determinant. As one can see, there is still 
a sign ambiguity remaining (originally we had two sign 
ambiguities) but that becomes irrelevant if we form the 
product of two pseudo-invariants. Indeed, the following 
quantity: 

^ _ Pf{<^(o,o,7r)}~^ det{L^(Q^Q^o)^(Q^o^^)}Pf{ (9(0^0^0)} 

\/det{[/(0,0,-7r)^(0,0,7r)} 
Pf{^(0,7r,7r)}~"^ det{?7(Q^^^0)^(Q^^^^)}Pf{^(Q^^^Q)} 

X , : 

vdet{C/(0,7r,-7r)^(0,7r,7r)} 



is a genuine invariant taking the quantized values ±1, 
which are independent of the branch of the square roots 
used in the calculation, as long as they are chosen consis- 
tently using the interpolating procedure described above. 
We can repeat the same construction for the pair of paths 
Vtt^o and Vn^TT and define the invariant: 

^ Pf{^(7r,0,7r)}~^ det{?7(^^Q^Q)^(^^Q^^)}Pf{6^(^^Q^Q)} 

\/det{t/(7r,0,-7r)^(7r,0,7r)} 
Pf{<^(7r,7r,7r) }~"^ ^^H^(7t,7t,0)^{7t,7t,7t) }Pf{^(7r,7r,0) } 

X z . 

V det{t/(7^^7^^_7^)^(7^^7^^7^) } 

The invariants Eq and are two of the four indepen- 
dent weak Z2 invariants. We can define two more weak 
invariants by pairing the paths in different ways, but that 
is not necessary because at this point we can define the 
strong Z2 invariant as: 

S = SoS, (25) 

If we count the strong Z2 invariant, then there are only 
3 independent weak Z2 invariants remaining. We will 
concentrate entirely on the strong invariant. 

We have already discussed the numerical aspects re- 
lated to computing the Z2 pseudo-invariants for each of 
the four paths of Eq. 20. There is another important nu- 
merical aspect about determining the correct branch of 
the square roots. One should note that computing the 
pseudo-invariants involves one dimensional calculations, 
in the sense that we only need to integrate along the 
kz direction and not on a surface as it is the case when 
applying, for example, the popular algorithm of Fukui 
et al from Ref. 33. However, we still have to perform 
the interpolation along the ky direction, so the calcula- 
tions become 2- dimensional. The key observation is that 
the number of ky points required by a successful inter- 
polation is usually an order of magnitude smaller than 
the number of kz points needed in the computation of 
the pseudo-invariants. This is because all we need is to 
determine how the determinants wind around the origin 
during the interpolation and to trace these paths one can 
indeed use a relatively small number of ky points. There- 
fore, our algorithm, while not strictly 1-dimensional, it 
can be regarded as quasi-one dimensional. 

To summarize, the application to the disordered sys- 
tem was possible because of the following numerical ad- 
vantages of the present algorithm: 

1. The algorithm is gauge independent. Finding a 
smooth gauge for a unit cell containing thousands 
of quantum states would have been practically im- 
possible. 

2. The quantization of the pseudo-invariants remain 
exact when the paths are discretized, allowing a 
drastic reduction of the number of the discretiza- 
tion points. 

3. The interpolation between the different time- 
reversal invariant paths can be accomplished with 
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FIG. 2. (a-c) The band structure of the model for a slab 
configuration < ns < 30, plotted as function of ki with /c2 
fixed at A:2 = 0. The hopping parameter t takes the values 
t = 14 meV in panel (a), t = 22.6 meV in panel (b) and 
t = 40 meV in panel (c). Panel (d) reports a calculation of 
the strong Z2 invariant as t was varied from 14 to 40 meV. 



a small number of k-points, transforming the algo- 
rithm into a quasi-one dimensional one. 



III. THE MODEL 

The model used in our numerical simulations is an ef- 
fective lattice Hamiltonian fitted to the topological ma- 
terial Bi2Se3. The starting point is a Hamiltonian Hq 
which, in the clean limit, can accurately describe the em- 
pirical energy band spectrum around the insulating gap. 
This Hq was used in the previous studies of disordered 
Bi2Se3 in Refs. 23 and 24. The Hq has inversion sym- 
metry and, since we want to exemplify the algorithms 
for systems without such symmetry, we will include an 
additional term in the Hamiltonian that strongly breaks 
the inversion symmetry. This term can be thought as the 
effect of a mechanical strain applied along the z axis. 

In the momentum space: 



Ho{k) = d^{k) 



( do{k) ds{k) 

d^{k) -do{k) d-{k) 

d^{k) do{k) 

\d^{k) -ds{k) 




where 



do{k) = e — 2t cos /c^, 
di{k) = —2X sin ki, i = 1,2,3 
d^{k) = 27 (3 - cos h) 



(27) 



d±(k) =di{k)±id2{k). 



(28) 



The added term that preserves the time-reversal symme- 
try but breaks the inversion symmetry is: 



Vt = R 



( 
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The following parameters will be fixed at these values 
throughout the paper: e = 134 meV, A = 30 meV, 7 = 16 
meV, R = lb meV. We will use t = 40 meV for the topo- 
logical insulator and t = 14 meV for the trivial insulator 
(the two values lead to comparable insulating gaps). The 
insulating gap in our study is larger than the empirical 
insulating gap of the Bi2Se3 material, and the reason we 
chose to proceed this way was to be able to better show- 
case the behavior of the strong Z2 invariant in the pres- 
ence of disorder (the insulating and the mobility gaps 
would have closed too fast if the gap was fixed at the 
empirical value). This modification does not break the 
bridge with the experimental reality because it is known 
that a mechanical strain may increase the insulating gap 
of the material. 

The real space representation of the translational in- 
variant Hamiltonian can be constructed on cubic lattice 
where each vertex n carries four quantum states |n, a, a). 
Here, a = ±1 (= isospin) labels the s or the p angular 
momentum character of a state the bands and <j = ±1 
the spin up and down configurations. On the Hilbert 
space spanned by |n,a,cr), we define <ii,j,/c, cr, a, and 

as the translation, spin, isospin and flipping operators 
as follows: 

dij^k\ni,n2,ns,a,(T) = \ni + i,n2 + j,n3 + k,a,a) 
a\n, a, a) = a\n, a, a), a\n, a, a) = a\n, a, a)r 
rcr\n,a,(j) = |n,a, — <j), ra\Ti,a,(j) = |n, — a,<j) 

(30) 

Then, the real space representation of Hq takes the form 
i^o = + 67 + A ^ sdo,s,ora^rcr 

s=±l 

-t{a + 7) X] (4,0,0 + f^o,s,o + ^^o,o,s) 

The term breaking the inversion symmetry takes the 
form: 



Vi 



a{a - l)(do,o,i^a " do.o.-O^o 



(32) 



As it is now well established, the topological properties 
of the clean model are revealed when restricting the total 
Hamiltonian Hq -\-Vx on di slab: < ns < N, where N is 
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taken large enough so that the tunneUng between the two 
surfaces of the slab is negligible. The slab Hamiltonian 
takes the form: 



-7 ^ (doos) + d'^i^i, /C2) + d^iki, k2)a (33) 

s=±l 

^di(ki)rara - id2{k2)ra^ra 

where d'^{ki^ /C2) = 27(3 — 2 cos(/ci)) and <io(/ci, /C2) = 
e — 2t^-y 2Cos(/c^). For such configuration, the parallel 
component to the surfaces of the momentum is conserved, 
so one can plot the energy spectrum as function of ki 
and k2- In Fig. 2(a-c) we show sections of such plots, by 
holding /c2 at /c2 = 0, for three different values of t. The 
dense band spectra seen in all three plots correspond to 
the bulk and one can see a bulk energy gap in panels (a) 
and (c). The bulk gap is closed in panel (b) and that 
marks the transition from the trivial to the topological 
insulator. Indeed, in panel (c) one can observe chiral 
bands connecting the valence and the conduction bands, 
and in panel (a) these bands are missing entirely. The 
chiral bands in panel (c), if plotted as function of ki and 
/c2, will give rise to a Dirac cone. The transition point 
between the phases is at t = 22.6 meV. 

A straightforward test of the algorithm described in 
the previous section consists of computing the strong Z2 
invariant for the clean system as function of parameter 

and comparing the output with the appearance or dis- 
appearance of the surface states in the slab calculations 
in Fig. 1. Fig. 1(d) reports these calculations and indeed 
both the Z2 invariant and the slab calculations predict a 
trivial insulator for t < 22.6 meV and a topological insu- 
lator for t above this value. The strong Z2 invariant was 
computed using the twisted boundary conditions on a 
4x4x4 lattice. The size of the lattice is irrelevant for the 
clean systems, and we just chose a convenient lattice size 
in order to test the twisted boundary conditions method. 
We used X number of kz points to compute the pseudo- 
invariants, in which case | det{/77r,-7r}| ^ ZZZ^ and Y 
number of ky points to perform the interpolations. 



IV. EFFECT OF DISORDER: LEVEL 
STATISTICS ANALYSIS 

The total Hamiltonian will include a random potential: 
H = Ho^Vx^V^, (34) 
where K; is a non-magnetic random potential: 

V^ = W ^ uJnoc\n,a,(7){n,a,(7\ (35) 

with ujjioc random entries uniformly distributed in the 
interval [— ^, 
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FIG. 3. Level statistics analysis for (left column) the topo- 
logical insulator t = 40 meV and (right column) the trivial 
insulator t = 14 meV. Each panel displays the variance of 
the level spacings ensembles as function of the energy where 
the level spacings were collected. The gray lines in each pan- 
els represent the integrated density of states (IDOS), which 
can be used to assess the evolution of the spectral gap, corre- 
sponding to the flat IDOS, and especially to determine when 
the gap is closing and becoming completely filled with local- 
ized states. The horizontal dash lines mark the value 0.104, 
the variance of the GSE ensemble. The vertical range in each 
panel goes from to 1. The shaded regions represent the 
emerging phase diagram of the topological model. 



For level statistics analysis, we diagonalized the dis- 
ordered Hamiltonian on a 14x14x14 lattice with peri- 
odic boundary conditions and for 500 random disorder 
configurations. We sampled the energy spectrum at 100 
equally spaced energy levels E. For each such we iden- 
tified, for each disorder configuration, the unique energy 
levels Ei and ^i+i satisfying: Ei < E < and we 
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recorded the level spacings: AE = — let- 

ting j take consecutive values between —5 and 5. Note 
that j indexes the levels and that each level is doubly 
degenerate. In this way, we have generated ensembles 
containing 5500 level spacings for each energy E. Fig. 3 
reports the variance (5^)/ (5)^ — 1 of these ensembles as 
function of energy E and disorder strength W. It also 
reports the integrated density of states (IDOS), which 
counts the number of eigenvalues below an energy E and 
normalizes this number by the dimension of the Hilbert 
space. When plotted as function of E^ the IDOS remains 
flat in the spectral gaps, so it is a useful and effective 
tool for identifying the spectral gaps in the energy spec- 
trum. We will be particularly interested to see when the 
insulating gap is closing. 

The level spacings follow a Poisson distribution when 
E is in the localized spectrum and the localization length 
is smaller than the size of the system. The Poisson dis- 
tribution has a variance equal to 1. In the spectral re- 
gions where the localization length exceeds the size of 
the system, the statistics of the level spacings coincides 
with that of a random Gaussian Symplectic Ensemble 
(GSE):4^'46 Pgse(5) = ^s^e-l^'\ The variance of 
this distribution is 0.104. One can study the trends as 
the system size is increased, and if the size of the system 
reached a limit where the variance is seen to stabilize 
(which we have verified that it does), the level spacing 
analysis can be used quite effectively to identify the re- 
gions of localized and de-localized spectrum. This will 
be done after we discuss the qualitative behavior of the 
energy spectrum in response to disorder. 

In 3 dimensions, extended states can exist even in triv- 
ial disordered models. The qualitative behavior of the 
spectrum in trivial models is as follows. Usually, the 
edges of the bands starts to localize the moment the dis- 
order is turned on (for systems displaying large variations 
in the density of states, additional patches of localized 
spectrum can occur deep inside the band). At moderate 
disorder, extended states still survive deep within the 
bands. As such, there are usually two mobility edges 
forming per band, flanking the region of extended states, 
and these mobility edges moves towards each other when 
the disorder is increased until they merge and disappear. 
At that point, all the states become localized, as it should 
be the case at large disorder. ^'^ 

In a topological model, the behavior of the spectrum 
is markedly different. The edges of the bands are still 
the first parts of the spectrum to become localized and 
the extended states are still located in the middle of the 
bands. But, if a band carries a nontrivial topological 
number that is robust against disorder, the two mobility 
edges flanking the extend states in a band cannot merge 
and disappear like in the trivial case because that will 
lead to a sudden change of the topological number car- 
ried by the bands, from a nontrivial to a trivial value. 
So what happens when increasing the disorder? The 
energy spectrum will eventually become entirely local- 
ized as the disorder is being steadily increased, and the 
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FIG. 4. Illustration of the phase diagram of the model as 
derived in Fig. 3, and the two paths used in the mapping of 
the Z2 invariant. 



only way this can happen is through a scenario where 
the bands carrying topological numbers collide with each 
other and in the process they neutralize their topologi- 
cal numbers. This leads to one of the hallmarks of the 
topological models where, when increasing the disorder, 
the spectral regions of extended states are seen to drift 
towards each other until they merge and disappear, usu- 
ally at very large disorder strengths. The levitation of 
the Chern-number carrying extended states in the Inte- 
ger Quantum Hall Effect is well known from the works 
of Halperin and Laughlin,^^'^^ and the pair annihilations 
of the topological states in lattice models of IQHE was 
discussed in Refs. 50 and 51. The levitation and pair an- 
nihilation picture was instrumental for the understanding 
of the global phase diagram of IQHE, and that will also 
be the case for our study. 

In Fig. 3 we report the variance of the level spacings 
ensembles collected at various energies and for increas- 
ing disorder strengths. We do not show here the actual 
histograms of the level spacings because of the large vol- 
ume of data already contained in this figure. However, 
the level statistics have been exhaustively researched for 
topological models, ^^'^^'^^'^^ and the correlation between 
the histograms and the value of the variance has been al- 
ready firmly established.^^ Panels (al)-(all) refer to the 
topological case where t = 40 meV and panels (bl)-(bll) 
refer to the trivial case where t = 14 meV. These two 
values were chosen such that the insulating energy gaps 
are practically the same (see Fig. 1). 

Examining the panels in Fig. 3, one can observe energy 
regions where the variance is large (and becomes unity at 
large disorder) but also energy regions where the variance 
remains pinned at the 0.104 value. These later regions 
will be identified with the spectral regions of extended 
states, while the former ones with the spectral regions 
of localized states. In panels (al)-(all) we can clearly 
see the two extended states regions drifting and merging 
with each other as the disorder is increased. The ex- 
tended states survive even at extreme values of disorder 
W = 1000 meV; this value is about twice the width of 
the entire clean energy spectrum. No such behavior is ob- 



9 



served for the trivial case in panels (bl)-(bll), where the 
valance band is seen to become entirely localized already 
at VF's as small as 200 meV, and the whole spectrum 
becomes localized before W reaches 700 meV. 

Based on the data presented in Fig. 3(al)-(all), we can 
draw the phase diagram of the topological model in the 
(VF, Ep) plane with quite accurate precision. It consists 
of a strong topological insulating phase surrounded by a 
metallic phase, which at its turn is surrounded by a triv- 
ial insulating phase (see Fig. 3). One could be inclined, 
by just looking at this 2D phase diagram, to call this 
later phase the Anderson insulating phase rather than 
the trivial insulating phase, but this is not correct be- 
cause if we consider a third dimension to the phase dia- 
gram along the parameter t, one will easily see that this 
phase is connected to the trivial insulating phase, such 
as for example the t = 14 meV and w = case shown in 
Fig. 2. As such, this phase should be called simply the 
trivial insulating phase. Now, examining the integrated 
density of states, we can see that the spectral gap is al- 
ready closed at 1^ = 300 meV but the topological phase 
extends beyond this W value. This phase diagram will 
be re-confirmed by a direct mapping of the Z2 invariant. 



V. MAPS OF THE Z2 INVARIANT 

The Z2 invariant will be mapped along the two paths 
shown in Fig. 4. Due to the extreme computational costs 
of such calculations, we had to settle for a somewhat 
smaller lattice size of 8x8x8 (but same as the largest 
lattice size used in Ref. 23). We have used 400 /c-points 
in the kz direction to compute the monodromies, and 
25 /c-points in the ky direction for the interpolation. As 
such, each Z2 invariant computation requires 20,000 ex- 
act diagonalizations of the disordered Hamiltonian. A 
number of 10 disorder configurations were considered for 
each point chosen along the paths shown in Fig. 4. 

There is one important numerical aspect that we must 
acknowledge, which is the computation of the pfaffian of 
the time-reversal operator at the time-reversal invariant 
A:-points. This became an issue for us because the dimen- 
sions of the matrices are very large. We have successfully 
used the fortran routine PfaffianH freely provided by the 
authors of Ref. 56, which computes the pfaffian of a gen- 
eral complex a skew- symmetric matrix using the House- 
holder transformations. 

Fig. 5 reports the map of the Z2 invariant along the 
path (1). The calculations were performed with a fixed 
number of occupied states rather than a fixed Fermi level. 
By doing so, we ensured that all the projectors in the 
monodromy formula Eq. 14 have the same dimension but 
in this case the Fermi level displays small fluctuations 
which disappear in the thermodynamic limit. The dimen- 
sion of the occupied states was slowly reduced from 1024 
(half- filled) to 50, as illustrated in Fig. 5, and for each 
dimension we have computed the average Fermi level, 
defined as half between the last occupied and lowest un- 



occupied states. The averaged Fermi levels are shown 
as vertical lines, over-imposed on the variance plot at 
W = 300 meV. As one can see, the Fermi levels sample 
the entire spectrum below the gap. 

We want to point out again that the spectral gap is al- 
ready closed at = 300 meV but, according to the level 
statistics analysis, there is still a mobility gap opened. 
We have verified this statement by direct check of the 
eigenvalue files. Also, the integrated density of states 
shows an inflection point rather than a plateau. When 
the Fermi level was inside this mobility gap, we found ab- 
solute no fluctuations in the Z2 invariant, which turned 
out to be —1 for all 10 random configurations. As the 
Fermi level is lowered, it enters the region of extended 
states and here we observed large fluctuations. As we 
already discussed, even in this regime the Z2 continue to 
take quantized ±1 but there is no way to tell which one 
will be so the output fluctuates between the two allowed 
values. This remains the case for as long as the Fermi 
level is in the region of extended states and, as soon as 
the Fermi level emerges back into the region of localized 
states, the Z2 invariant is seen to take a non- fluctuating 
quantized value of +1. 

A similar behavior is observed when considering the 
path (2) of Fig. 4 for which the calculations are re- 
ported in Fig. 6. While increasing the disorder strength, 
the Z2 invariant is seen to take the quantized and non- 
fiuctuating value of —1 until the mobility gap closes. 
From there on, the values fluctuate between ±1, and the 
Z2 stabilizes once again when the path enters the trivial 
insulating state where it assumes the value -hi. 



VI. CONCLUSIONS 

In conclusion, a previously introduced gauge- 
independent formulation of the strong Z2 invariant was 
found to bring significant numerical advantages, allow- 
ing direct computations of the invariant for large super- 
cells with twisted boundary conditions. The resulting 
numerical algorithm was applied to a disordered model 
of Bi2Se3 topological material and maps of the strong 
invariant were given as function of either Fermi level or 
disorder strength. The behavior of the strong Z2 invari- 
ant seen in our numerical calculations is exactly what one 
will expect if this invariant was indeed robust to disorder. 
Specifically, we observed the strong Z2 invariant taking 
quantized and non-fluctuating values whenever the Fermi 
level was in an energy region of localized states, and fluc- 
tuating values (between the only two possible values of 
±1) whenever the Fermi level was in an energy region of 
delocalized states. The fact that our numerical maps of 
the strong Z2 invariant were in good agreement with the 
phase diagram constructed from the level statistics anal- 
ysis leaves very little doubt that the strong topological 
phase survives beyond the point where the spectral gap 
closes, and that it extends all the way to the point where 
the mobility gap closes. 
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FIG. 5. The upper panels show the invariant computed along the path (1) of Fig. 4, on a 8x8x8 unit cells lattice via 
twisted boundary conditions. The dimension of the occupied space was slowly reduced from 1024 to 124, as indicated in the 
figure. Each Zi the calculation was repeated for 10 random disorder configurations and the output is shown by the full dots, 
exactly how it occur in the actual calculation. The percentages of the = ±1 occurrences is displayed in each panel. The 
lower panel shows the variance of the level spacings for W — 300 meV, and the averaged Fermi levels (see the vertical lines) 
corresponding to each calculation. 



Our algorithm, combined with accurate tight binding 
models that can be developed for any material via ei- 
ther first principles calculations or by simple empirical 
means, can provide accurate quantitative simulations 
of the real experimental samples. We want to point out 
that, recently, the twisted boundary conditions were suc- 
cessfully used to compute the Chern invariant of an inter- 
acting 2-dimensional fractional Chern insulator. Since 
the algorithm for computing the Z2 invariants is less de- 
manding than the algorithm for the Chern invariant, we 
have high hopes that we will soon be able to map the Z^ 
invariants in the presence of electron interaction for accu- 
rate complex models of topological materials. Both dis- 
order and electron interactions are expected to strongly 
influence the phase diagram of a topological material. 
The topological/non-topological state of a sample can be 
probed by looking at the extended/localized character 
of the surface states via transport simulations on quasi- 
one dimensional bars (like it was done in Ref. 24), but 
here one has to be careful with the boundary conditions 
and, in addition, the extended surface states can exist in 
non-topological samples. Therefore, we believe that the 
direct computation of the strong Z2 invariant will be a 



valuable complement to these aforementioned methods. 

The data generated by our study can be of interest for 
experimentalists. So far, all topological materials fabri- 
cated in the labs display metallic bulk properties, a fea- 
ture that was attributed to the imperfections of the ma- 
terials. Our study revealed the interesting fact that, due 
to the very topological nature of the materials, the disor- 
der pulls the valance and the conduction mobility edges 
closer to each other. In fact, within our tight-binding 
model for Be2Se3, we saw a rapid reduction of the mobil- 
ity gap with disorder and the closer of the mobility gap 
when the disorder strength reached about 350 meV. This 
suggests that the topological materials have to be much 
"cleaner" than their trivial counterparts in order to see 
an insulating bulk phase. 
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